What can one learn about Self-Organized Criticality from Dynamical Systems theory ? 
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We develop a dynamical system approach for the Zhang's model of Self-Organized Criticality, 
for which the dynamics can be described either in terms of Iterated Function Systems, or as a 
piecewise hyperbolic dynamical system of skew-product type. In this setting we describe the SOC 
attractor, and discuss its fractal structure. We show how the Lyapunov exponents, the Hausdorff 
dimensions, and the system size are related to the probability distribution of the avalanche size, via 
the Ledrappier- Young formula Q . 
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I. INTRODUCTION. 

Within the last 10 years the notion of Self-Organized Criticality (SOC) became a new paradigm for the explanation 
of a huge variety of phenomena in nature and social sciences. It's origin lies in the attemp to explain the widespread 
appearence of power-law like statistics for characteristic events in a multitude of examples like the distribution of the 
size of earthquakes, 1/f-noise, amplitudes of solar flares, species extinction .... to name only a very few cases ]l|-p|,^9[. 
g ; As a result, an important literature in physics has been devoted to the study of systems exhibiting SOC. 

The complexity of the dynamics in the above mentionncd systems is mainly due to the presence of long-range spatial 
and time correlations, leading to non trivial effects like anomalous diffusion. At stationarity, the average incoming 
■ flux of external perturbations is simply compensated by the average outgoing flux that can leave the system at the 
tyQ " boundary, or by dissipation in the bulk. Therefore, there is a constant flux through the system, leading to a non- 
| equilibrium situation. What is remarkable in this stationary state, refered to as the SOC state, is that the distribution 
of avalanches appears to follow a power law, namely there is scale invariance reminiscent of thermodynamic systems 
at the critical point. This is certainly one central reason why SOC has attracted the physicist community: these 
systems (apparently) reach spontaneously a critical state without any fine tuning of some control parameter. 

Several models have been proposed to mimic these mechanisms like the sandpile model the abelian sandpile 
|l2| or the continuous energy model |3Q]. Numerical simulations on one hand, and theoretical approaches on the other 
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O ■ hand have lead to a good description of SOC, in particular with respect to the computation of critical exponents that 



are believed to characterize the universality class the model belongs to, as they do in second order phase transitions. 

However, to our knowledge, no serious attemp has been made to study SOC from a dynamical system point of view 
(except |^,|ll|]). It is however a natural approach to try to access the macroscopic behaviour of large sized systems 
from the microscopic dynamical evolution. The macroscopic behaviour at stationarity is characterized by a probability 
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measure one has to extract from the microscopic evolution. One is seeking a "good" measure from a physical point of 
view, namely a Sinai-Bowen-Ruelle measure (SBR) : in the SOC model we discuss later this measure maximizes the 
entropy. 

In this paper we develop a dynamical system description for a certain class of SOC models (like the Zhang's 
model pot), for which the whole SOC dynamics can either be described in terms of Iterated Function Systems, or as a 
piecewise hyperbolic dynamical system of skew-product type where one coordinate encodes the sequence of activations. 
Several deep results from the theory of hyperbolic dynamical systems can then be used, having interesting implications 
on the SOC dynamics, provided one makes some natural assumption (like ergodicity) which will be partially justified 
in this paper. 

With this approach we give a precise definition of the SOC attractor discussed by some people JlL§]. We show 
that it has a fractal structure for low values of the critical energy. The main objects for which our point of view is 
appropriate is certainly the structure of the asymptotic energy distribution or, in other words, the structure of the 
natural invariant measure. We show in particular how the Lyapunov exponents, the geometric structure of the support 
of the invariant measure (Hausdorff dimensions) , and the system size are related to the probability distribution of the 
avalanche size, via the Ledrappier- Young formula po| . 

II. THE DYNAMICAL STRUCTURE OF THE ZHANG MODEL. 
A. Description of the model. 

In this paper we deal with the Zhang's model on a d dimensional, connected subgraph A C 7L a , with nearest 
neighbours edges, though the formalism we develop holds for more general graphs. Let 9 A be the boundary of A, 
namely the set of points in 7L d / A at distance 1 from A and let N the cardinality of A. Each site i € A is characterized 
by its "energy" Xi, which is a non-negative real number. The "state" of the network is completely defined by the 
configuration of energies X = {A^} ieA . Let E c be a real, positive number, called the critical energy, and M = [0, E C [ N . 
Let <ii ) 2(X, Y) be the L\ (resp. L^) distance on M. A configuration X is "stable" iff X £ M. and "unstable" or 
" overcritical" otherwise. If X is stable then we choose a site i at random with probability -4, and add to it energy 
SX. As far as the physically relevant parameter is the local rigidity [^3|, one can investigate the cases where E c 
varies, and where 6X is a constant. We will therefore assume that SX = 1 without loss of generality. If a site i is 
overcritical (Xi > E c ), it loses a part of its energy in equal parts to its 2d neighbours. Namely, we fix a parameter 
e € [0, 1[ such that, after relaxation of the site i, the remaining energy of i is eXi, while the 2d neighbours receive the 
energy ^ 1 ~' 5 t ] X ' . Note that in the original Zhang's model pp| , e was taken to be zero. We define here a straightforward 
extension. Note however that in this paper e will be considered as a small parameter compared to E c . 

If several nodes are simulaneously overcritical, the local distribution rules are additively superposed, i.e. the 
time evolution of the system is synchronous. The sites of 9 A have always zero energy (dissipation at the bound- 
aries). The succession of updating leading an unstable configuration to a stable one is called an avalanche. 
Because of the dissipation at the boundaries, all avalanches are finite. The structure of an avalanche can be 
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encoded by the sequence of overcritical sites A = {Ai} Q<i where Aq = {a}, the activated site, and Aj = 
{j € A|Xj > i? c in the ith step of avalanche} , i > 0. 

The addition of energy is adiabatic. When an avalanche occurs, one waits until it stops before adding a new energy 
quantum. Further activations eventually generate a new avalanche, but, because of the adiabatic rule, each new 
avalanche starts from only one overcritical site. 

Since the avalanche after activation of site a maps overcritical to stable configurations one can view this process as 
a mapping from M. — ► M. where one includes the process of activation of site a. We hence associate a map T a with 
the activation at vertex a. This map usually has singularities and therefore different domains of continuity denoted 
below by M„ where k runs through a finite set depending on a. Call T% = T a \j^h. The main object of this paper is 
the study of the properties of the family of mappings {T^} and to link these properties to the asymptotic behaviour. 

B. Piecewise afflne mappings. 

1. Structure of the piecewise affine mappings. 

One can easily write the conditions on the stable energy configurations insuring that the avalanche A = {Ai} Q<i 
occurs. This defines a convex domain^] in A4 . The A4^s are the domains of continuity of T a and they constitute, 
for each a, a partition of M.. There is therefore a one to one correspondence between an avalanche and a map T*. 
The energy distributions rules of Zhang's model implies that: 

T*.X = L k a .(X + e a ) = L k a .X + L k a .e a , XgMJ. (1) 

where the linear mapping L k characterizes the redistribution of energies on each sites after the avalanche. The 
column i of s contains the ratios of energy given by the site j to the other sites after the corresponding avalanche. 
Alternatively, the entries of the row i correspond to the energy received by i from the others sites (i included). e a 
being the canonical basis vector of IR^ in the direction corresponding to the activation at site a, the constant vector 
L k .e a corresponds to the redistribution of the additional energy SX = 1 on each site, after the avalanche. In the case 
where no relaxation occurs the corresponding map is just a shift along the a axis. A way to build L\ is to construct 
it step by step, by a left product of elementary matrices giving the redistribution of energy from one step in the 
avalanche to the successive step. The composition of these matrices is determined by the avalanche profile. 

Let S% = dAia- Then S a = \J k is the the set of singularities for the transformation T a . The sets S% are unions 
of segments of hyperplanes in IR^ . 



x For example, the domain of energy such that a site a do not relax after activation is delimited by the hyperplane E a = E c — 1, 
and the boundaries of A4. 



3 



2. Projection effect. 



The original Zhang's model contains a pathology. Due to the reset to zero of an overcritical site after relaxation 
(e = 0), linear dependences among the sites (more precisely, direction of IR^ associated to these sites) are created 
along the avalanche. This implies the existence of a non trivial kernel. Thus each L^'s is a projection onto a subspace 
of R , whose dimension increases with the number of involved sites in the avalanches. The one step matrices have 
a number of zero eigenvalues given by the number of sites set to zero at the corresponding time step. Multiplying 
these matrices gives raise to the kernel of L\. We denote the subspace generated by zero eigenvalues £°(a, k). Note 
that, in general, KerL k a C £°(a, k) . Clearly, the existence of a nontrivial kernel is the source of several mathematical 
complications when studying the dynamics of Zhang's model. It is a very particular feature of the e = model. It 
makes however the global geometry of the attractor quite interesting (see Fig. [I]). 



3. Contraction 



Each mapping L k a has only eigenvalues of modulus lower or equal than 1. Indeed, by definition, the sites of the 
boundaries of the avalanche receive energy without relaxation. This implies that, by eventually permuting the basis 
vectors, the linear map L can be written as: 

I * ... * 

L= (2) 

* ... * I 

where I is the identity matrix and where the *'s can be zero or not. They correspond to the fraction of energy received 
by the sites which have not relaxed. Therefore, the vectors corresponding to sites not relaxing are eigenvectors of L 
with eigenvalue one. We denote the corresponding (neutral) subspace by £ n (a,k). 

On the other hand, the inner block in (2) corresponds to the sites which have relaxed. The energy conservation 
implies that the sum on each column of the block is strictly lower than one (some part of the energy has gone outside 
the block, to the sites on the boundary of the avalanche). By usual arguments on positive matrices it follows that the 
eigenvalues are strictly lower than one in the block Jl^| . Note that, for e = 0, this block contains also the subspace 
£°(a, k). Therefore, the subspace of relaxing sites is decomposed into two subspaces: £°{a,k) and £~{a,k), where 
£~(a,k) denotes the subspace associated to the eigenvalues < |A,| < 1. 

Hence to each mapping L\ we associate the following decomposition: 



£~ (a, Jfe) © £ n (a, k) © £°{a, k) = IR 



N 



(3) 



C. Composed mapping. 

1. Composition of affine mappings. Extended dynamical system. 



The activation dynamics can be represented by the left Bernoulli shift a over , the set of right infinite sequence 
a = {<zi, . . . , dfc, . . .} , afc £ A, where era = 02(13 . . .. Namely, a n is the n th activated site in the activation sequence a. 
We denote by [a] the set of sequences whose first digit is a. 
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The combined effect of the activation and relaxation process is than described by a dynamical system of skew- 
product type T : — ► £1 such that: 

T(X) d ^>a,T Ql (X)) ; X%\a, X ) 

where f2 = x A4 is called the extended phase space. 

Let DT-^ be the tangent map of T at X. (When speaking about differentials of T we usually think of Sa represented 
by a smooth system z — * |A| • zmod 1). 

The singularity set of T is : 

S = |J [a] x S a (4) 

a€A 

We define a distance on f2 by g?q(X, Y) = d E + (a, a') + oIm (X, Y), where X = (a, X) , Y = (a', Y). We denote the 
two projections on the first and second coordinate by tt u (a, X) = a, and ir s (a, X) = X. The superscript it, s means re- 
spectively unstable and stable and will be explained below. We have a natural partition of il, V — \V\ — [a] x .M*}. 

Note that V is a generating partition for (T, f2) in the topological sense, that is, the diameter of the elements of 
Y _ f-ij) g 0es t zer0 _ 



2. Kernel of the infinite product map. 

For e = the kernel of the map T* can increase with i, projecting IR^ onto spaces of lower and lower dimensions. 
Therefore, after a certain, finite time, n(X), X is projected onto the effective lower dimensional subspacc ^| [J 

£ s (±) = {v g {0} x R N ; Vi > 0, ||2W£.v|| > o} (5) 

It is somehow the reference space with respect to X, because, asymptotically, the dynamics of vectors in IR W under 
DT-^ reduces to the dynamics of vectors initially in £ S (X). We get therefore a splitting of the projection on Ai of the 
tangent space at X as: 

IR^ =£ s (X)e/C(X) (6) 

where /C(X) is the kernel of the product map. This splitting will be refined further below by using the Oseledec space 
decomposition. 



2 One has to check that this space is not finally reduced to {0} ! However, this would imply that all vectors in are 
asymptotically mapped to 0. This is excluded since the avalanche matrices have non negative entries. 

3 Note that the dimension of this space depends a priori on X. However, if ergodicity holds, this dimension is constant for 
(almost-every) initial condition. 
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3. Local contraction. 



One easy shows that for any finite connected A and arbitrary activation sequence a all sites become overcritical 
infinitely often. Assume the opposite. Then, there exists a site which is overcritical only a finite number of times along 
the infinite sequence a but has a neighbour site which relaxes infinitely often. The energy coming from the overcritical 
neighbour site is larger than ^—£f^ by definition. This implies that all neighbours relax also an infinite number of 
times during the whole sequence. Hence we get a contradiction. It follows that there exists a time r = t(A, E c , e) < oo 
such that, VX, after at most r time steps each site has been at least once overcritical. By looking at the product map 
DT^ this implies that all eigenvalues are different from one. This is straightforward since the sum of entries on each 
column of the composed map on A4 is strictly lower than one (and is bounded away from 1). Therefore there is a 
positive constant C = C(A,E c ,e) s.t.: 

h'DTTW, = sup \\ifDTLV\\t < C < 1 . (7) 
\Wh=i 

This implies that the map T T acts as local contraction in all directions in the space A4, along the trajectory of any 
point X . This has in particular the following consequence. The distance of two points X, Y whose trajectory belong 
to the same domain of continuity eventually goes to zero if the trajectories lie in the same domains of continuity along 
the whole activation sequence. 

4- Hyperbolic structure and Lyapunov exponents. 

Assume that almost every point is regular, namely the map T is differentiable along all points of the trajectory 
(note that as long we work in the tangent spaces this assumption is not necessary since the involved mapping are 
all well defined at dPi. Only for the construction of the local induced stable manifolds in M one has to take care of 
regularity). Then, one can decompose the (effective) tangent space at a.e. point X £ f2 into a contracting subspace 
£ S (X) and an expanding one £"(X), s.t. : 

1. VX, £ s {±) © £ u (±) = R N+l-dim(IC(±)) _ 

2. VX, 3A < l,r < oo, s.t. \\DT T \\ £S{±) < A. Furthermore \\DT\\ £u{±) = N = |A|. 

3. T{£ s {±)) = £ S (T(X)); T(£ u (±)) = £ U (T(±)) 

Furthermore £ S (X) can be decomposed into a sequence of subspaces Q,|^,|9| 0: 

£ S (X) =£ 1 (X) D£ 2 (X) d ... D£i(X) (8) 
such that if v G £j(X) \ £; + i(X) the average contraction of v is given by the Lyapunov exponent : 



4 In fact, in the e = case we have still to assume that the angle between £ S (X) and /C(X) is bounded away from zero, a.s. 
because, otherwise, Lyapunox exponents might not exist 
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Ai(X) = lim -log\\DT£(y)\\ 3 (9) 

From property (^) there are no zero Lyapunov exponents (note however that some exponents go to zero as E c tends 
to infinity). 

If the dynamics is ergodic the Lyapunov exponents are almost-surely constants and the same holds for dim£i. 
Corresponding to the shift action there is a positive Lypunov exponent, which is trivially log(N). The Lyapunov 
exponents are directly related to the geometrical structure of the support of the invariant measure. In the Zhang's 
model the negative exponents > Ai > A2 > . . . > Xn are physically related to the transport of energy and the 
dissipation rate at the boundary ||. In particular, we show below that there is a natural relation linking these 
exponents to the avalanche size distribution. The positive Lyapunov exponent Ao = logN corresponds to the entropy 
production coming from the activation dynamics. 

The average exponential volume contraction rate on £ S (X) is given by the sum: 

N 

E A * ( 10 ) 

z=l 

while the average exponential variation rate of the volume in the extended phase space is logN + y'., A^. 
For regular X let 

W £ S (X) = {Y : d(±, Y) < e Vt > 0, T l X , 

and T 4 Y are in the same partition element of vX 

be the e local stable manifold. Clearly one has on W| uniform exponential contraction. 
The global stable manifold W S (X) is obtained by 

\Jt-*(w s a±) . 

i>0 

Finally let W; c(X) be the largest connected component of W s containing X. Since the system is of skew product 
type one has a trivial unstable manifold being in the case of representing the shift as z — > z ■ |A|modl the whole 
interval [0,1]. Note that W e s (X) may not exist if 3{ni} s.t. d(T ni (X),5) < e~™ iC where C > is some constant 
larger than — Ai. The set of points with this property has measure zero unless the invariant measure concentrates on 
S. This aspect will be described in more detail in ||. We make the following conjecture: 

Conjecture 1 There exists a E C (N), such that, for Lebesgue almost- every E c < E C (N) there exists an n(E c ,N) and 
a v such that Vt > n(E c , N): 

d(T'(n),5) > v > 

This implies that after a finite time the dynamics stays away from the singularity set. This assumption is sufficient 
for the existence of local stable manifolds, but it will have several other important implications. We expect Conj. II] to 
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be true for E c sufficently small since for E c <C 1 the contraction dominates the expansion in the extended phase space. 
In this case the invariant set has the structure of a totally disconnected Cantorset with large gaps. Furthermore the 
local structure of this invariant set is constant for open sets of E c values (see section 6), but the singularity set is 
varying continuously in E c (except for a countable set), one can find open domains of E c values where S stays away 
from the invariant set. 

The singularity set has nevertheless the following effect on the dynamics. Take an 77-ball of initial conditions (in 
M.), and fix an activation sequence. For rj large enough the image of the ball under some iterate of T is cut by the 
singularity set. This means that the points separated by the singularity set will not evolve under the same sequence of 
mappings. For large size systems this can cause on M. a kind of expansion effect (for fixed typical activation sequence) 
on a mesoscopic scale. 

5. Symbolic dynamics. 

Symbolic dynamics is a very useful tool for the investigation of the orbit structure of dynamical systems. To this 
aim, one fixes a partition V of the phase space and associates to each point the sequence of partition elements visited 
by the orbit of a point. To make symbolic dynamics useful one wants this correspondence essentially to be unique 
(that is up to sets of measure zero) . Furthermore to handle the symbolic dynamics it is essential to have an explicit 
characterization of the legal (that is by orbits generated) set of symbolic sequences. The perhaps most prominent 
example of such an explicit description are symbolic systems defined by a Markov transition graph called Subshift 
of Finite Type (SFT). A classical result in hyperbolic dynamics says that uniform hyperbolic systems are always 
conjugated to SFT p0| , p5| , p7[ . The specific partitions giving rise to such coding are called Markov Partitions. 

In the Zhang's model one can encode the possible transitions between avalanches in a transition graph with respect 
to the canonical partition V = {P„ }. Namely, we draw an arrow from to V\ if and only if T(V^) (~l V\ 7^ 0. We 
denote by the set of admissible infinite sequences w.r. to the partition V . . 

Clearly, points on W;* C (X) form an equivalence class for the symbolic coding induced by V . Note that the transition 
graph is a priori Markov only for special choices of E c . 

When the invariant set is bounded away from the singularity set one can refine the partition {'P*} to make it 
Markov (this is certainly not a necessary assumption to get a Markov transition graph). Namely, there is an m s.t. 
the partition \/™ =1 T~' l (T' m V) is a Markov partition. We label the affine mappings corresponding to the Markov 
partition elements by Fi. Note that several Fi can usually correspond to the same map T*. 

Let us give an example. In the case E c s [1,2], e — 0, in one dimension, the piecewise continuous mapping applied 
is uniquely determined by the position of the zero site. Indeed, after a sufficiently long sequence all sites have energy 
Xi > ^f- but eventually one with a zero value. Since a site with value zero is the only possible stopping site for an 
avalanche besides the boundary the avalanche is uniquely determined by the position of the activated site and of the 
zero. This case is however the simplest, because there is no need to cut further the V^'s in order to get a SFT. For 
E c > 2 things are more complicated, due to the presence of sites with integer values 1 . . . [E c — I] which may stop an 
avalanche, according to the amount of energy they receive. 
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As already said we expect the above mentioned property of being disjoint from the singularities to hold for a.e. E c 
value less than some E C (N). This implies that the system is a SFT. For e = 0, d — 1 there is another dense set of E c 
values for which one can show that the system is a SFT. 

Proposition 1 For e = 0, d = 1 and E c — n/(2d) p for any n,p 6 IN* the system is conjugate to a SFT 

Note that the elementary operations on each avalanche and each node i are of the form Xi —* Xj + Y^j ^ for some j 
and a check wether Xi is larger or less than E c . If E c is of the above form it is a finite digit number (eventually zero) 
in base 2. It follows that for each point X £ M. one has to know only a uniformely bounded, finite number of digits 
in base 2 to decide in which set M.\ X is. The same holds for the legal transitions between avalanches domains, 
that is, there is a finite number of forbidden strings in base 2' A ' coding the whole system, hence it is a SFT. □ 



6. Macroscopic state and SBR measures. 

The addition of energy on one hand, and the dissipation of exceeding energy at the boundaries, on the other hand, 
drives gradually the system towards a stationary state where there is a constant energy flux through the system. 
As far as our representation accounts for activation dynamics on one hand and transport-dissipation (avalanche) on 
the other hand, the full informations about the macroscopic behaviour of the system at stationarity is contained in 
the invariant measures^ of our dynamical system. Since Q has a product structure one has canonical measures /x M 
(induced measure on the unstable direction) and /i s (induced measure on A4). For simplicity we will assume that 
fi is a Bernoulli measure, namely that the sites are chosen independently with fixed rates. Once we have fixed the 
distribution of activation, we are interested on the possible fi s measures. Of special physical importance are the 
measures obtained by iterating the Lebesgue measure \il onM, that is lim^oo i J27=o T % {[l u X We call this 
measure conditional SRB with respect to fi u . 

It is common in the SOC litterature to assume ergodicity. In our setting the physically relevant ergodic property 
is equivalent to the following conjecture. 

Conjecture 2 For any E c ,A,e, and given fi w the corresponding conditional SBR measure is unique. 

This implies in particular the almost-sure equality between the ensemble average and the time average for typical 
energy configurations. We give some arguments to support this assumption at least for certain E c values. If the 
probability of activation of any site in non zero then there exists a periodic point X with period p such that 



/' 

where 



>0 (11) 



7r"(X), 7r n (T(X)), . . . , 7r"(T p_1 (X)) is a cylinder set. This is the set of infinite sequences in which coincide 
with the activation sequence of X on the p first symbols. Assume that the periodic orbit admits a stable manifold 
such that 



5 Namely )i(T *(£>) = fi(B) where B is a measurable set in fl. 
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|JW s (T l (X)) = M (12) 

i 

If the periodic point does not lie on S one can take a small neighbourhood ZY e (X). Otherwise, as far as the singularity 
set is moving with E c while the limit cycle does not change on open domains of E c (the maps remain the same) one 



can change E c by an arbitrary small value in order to make the limit cycle disjoint from S. Due to (11), a generic 
sequence a admits arbitrary long segments with repeated words 7r u (X), 7r"(T(X)), . . . , 7r"(T p_1 (X). Therefore, from 
(|l2|), almost every points visits W e (X) infinitely often for any e > sufficiently small. By the hyperbolic structure 
and some moderate assumptions on the distribution of the size of Wf oc (X) one can then form a Hopf chairf^] between 
iterates of a.e. points Z, Y when they visit £Y e (X). By standard arguments from ergodic theory concerning the equality 
of forward and backward averages one can then prove that a.e. pair of points on the invariant set belongs to the same 
ergodic component. In general, it does not seem easy to give explicit examples of sequences of avalanches satisfying 
the above conditions. The crucial point here is to show (12). But perhaps it should be possible to weaken the above 
conditions substantially. 

For d = 1,E C > 1, one can check it by using the following argument. For E c > 1, starting from any stable 
configuration, one can add energy to the low energy sites < E c — 1) in order to get a configuration where all sites 
have energy E c — 1 < E{ < E c . Activating any site in this configuration generates a unique "maximal" avalanche 
where all sites become overcritical. This avalanche is recurrent [] and there exists a periodic orbit satisfying (12). 
For d > 1 the number of reflexions of the front on the boundaries can vary and there are several types of " maximal" 
avalanches which makes the argument break. One can however still apply it on a diamond shaped lattice with with 
L odd (N = L 2 ) , because, by activating periodically in the middle site, one has essentially the same situation as for 
the one dimensional chain. For E c small, especially for E c < 1 the avalanche patterns are much more complicated 
and the above argument breaks down. 

Note that in the case where we have a Markov graph, the ergodic property can in principle be directly checked on 
the Markov transition graph defined in the previous section. Namely, if the Markov transition graph is asymptoticaly 
irreducible and aperiodic, (unique) ergodicity follows from usual results on Markov chains. 

We proceed in discussing some aspects of the dependence of \i on E c for fixed A. 
Proposition 2 fi s is singular for all E c sufficiently small. 

Proof This follows easily since for E c << 1 one can make the L\ norm of all avalanche map arbitrary small since the 
avalanche has to "reflect" many times on each boundary node, hence every node has contributed to the dissipation. 
Since the expansion is constant it follows that det(DT) < 1 hence all measures are singular. □ 

Proposition 3 fi s is atomic for the chain and E c £ [jzr, jzr]- 



6 A path made of pieces of local stable and unstable manifolds. 

7 This argument has been already used by other authors like Dhar p3|, and Speer [|8j for the Dhar model. 
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This is proved in section III. A. 3. Furthermore, we conjecture the following: 

Conjecture 3 The Hausdorff dimension of fj, s is piecewise continuous and monotonously increasing on the domains 
of continuity for E c « 1 . 

This is supported by the following argument. On open intervals Ii of E c the structure of the mappings T* does not 
change but the domains of continuity hA\ do. Furthermore for E c decreasing the probabilities for avalanches with 
higher contraction should increase which should force the Hausdorff dimension to increase monotonously with E c on 
each It. 

We now discuss the connection between the invariant measure and the SOC state. It is possible to extract from /x the 
probability distribution of all observables usually considered in the study of SOC. The traditionaly used observables 
are: the duration t (number of iteration steps inside one avalanche); the size s (total number of relaxing sites counted 
with multiplicity), and the area a (number of distinct relaxing sites). Fix now an observable, say s. Let K, s be the set 
of mappings T* with avalanche size s and let Q s be the union of it's domains hA k a . Let Pn(s) be the probability to 
have an avalanche of size s for a lattice of size N in the staionary limit. One has clearly: 

P N (s) = fi s (Q s ) (13) 



7. Ledrappier- Young Formula. 

This formula plays a key role in relating the probability of avalanche size to the average contraction rate (sum of 
Lyapunov exponents). It establishes a kind of conservation law relating the Lyapunov exponents, some version of 
Haussdorf dimension and the Kolmogorov- Sinai entropy. 

One can refine the foliation into a stable and unstable manifold by splitting the manifolds into sub-manifolds W? (X) 
(resp W l Il (X)) such that the contraction (resp. the expansion) on Wf(X) (resp. W l Il (X) ) is governed by the Lyapunov 
exponent A^. Let Si be the local Haussdorf dimension of the measure fj, projected on W[(X.) (where r stands for s, u), 
namely: 

logfi(Bi(±,e)) 

b l = hm e ^ (14) 

loge 

where £>i(X,e)) is an e-ball around X in yV, r (X). Then ui = Si — Si+i,i = 1 . . . N — 1 is the transverse dimension 
of the measure fi on Wf(X) \ W l r +1 (X). It is constant for fi almost-every X if fi is ergodic. The unstable foliation 
being one dimensional in our context, the Haussdorf dimension of the measure yV"(X) is Sq. It is equal to one for the 
uniform activation measure. 

Let hfj, be the Kolmogorov-Sinai entropy of fi and the positive Lyapunov exponents. The Ledrappier- Young 
formula is (2(]] (for ergodic measures): 
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where the sum is taken over the positive Lyapunov exponents. It expresses in particular that without any absolute 
continuity of fi, any equation relating entropy and positive Lyapunov exponents must involve some notion of fractional 
dimension. In our case, it reduces to : 

I hl = log N S (16) 
From now on we will assume that \i u — \il (uniform activation). In this case Sq = 1. 

When the dynamics is invertible, this formula, applied to the inverted system, gives the following equality in the 
Zhang's model: 

N 

Ai (Ti = -logN (17) 

i=l 

where the sum is now taken over the negative Lyapunov exponents. 

However, one has to assume that the dynamics is (/i almost-surely) invertible. That means physically that, at 
stationarity, the probability that two avalanches, starting from two different configurations, end on the same config- 
uration of energies is zero. Like conjecture [l] we expect this property to hold only for small E c values, where the 
invariant set is a Cantor set but to fail generically for large E c values. 

We have the following conjecture: 

Conjecture 4 For E c sufficiently small, there exists a n(E c , N, d) such that Vt > n(E c , N): 

M (T*(^)nT t (Pj)) = 0, Vif^Pj 

Note that one can weaken this assumption by requiring that there are, on the attractor, less than N preimages and 
still get a nontrivial relation to Lyapunov exponents. One can still write down a Ledrappier- Young formula for non 
invertible systems by making the system invertible p6| by coding the backward iteration tree in the same way as we 
did with the activation sequences, hence introducing an additional variable on which the forward dynamics contracts. 
Let Jjv(X) be the number of preimages of X and Jjv = / Jj\r(X)d/i(X) the averaged number then: 

N 

— A-jO, = logN — logJN (18) 



III. DYNAMICS AND SOC. 
A. The Zhang's model as an iterated function system. 

If the system is conjugate to a subshift of finite type, the dynamics of the Zhang's model is essentially equivalent 
to a graph probabilistic Iterated Function System (IFS) ||[l6]], namely, a set of quasi-contractions i*j randomly 
composed along a Markov graph admitting a unique invariant measure fi* . Note that IFS are usually defined for true 
contractions, however, in our case, any finite composition along the graph is a contraction. In this case, the classical 
theory of graph directed Iterated Functions Systems applies and allows one to obtain interesting results with respect 
to the geometrical structure of the invariant set. 
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1. The Zhang's model attractor. 



The IFS determines a unique non-empty compact set A, called the attractor of the IFS, satisfying: 

dcf 



A = T{A) = {jF l {A) 



(19) 



This set is usually a fractal. 

Let H(A'l) be the set of compact subsets in M.. Define a distance on IH(A^), called the Haussdorf metric by: 

S(A, B) = sup {d(a, B),d(b, A),aeA,be B} 



(20) 



where A, B are non empty closed bounded subsets of A4, and d(x, A) = inf {d(x, a), a 6 A}. A is an attractor of the 
IFS in sense that it satisfies the following property Q : 

VB g W{M),F*{B) -» A 

in the Haussdorff metric when n — ► oo. Furthermore, if £> € H(A^) is such, that for all i, Fi(B) C £> then : 

oo 

^ = n f n {B) 

n=0 

Therefore, the asymptotics dynamics of the Zhang's model lives onto an attractor, further on denoted by A, whose 
fractal geometry is linked to the critical behaviour at stationarity. Note however, that, despite one might expect from 
the presence of dissipation the existence of an attractor with a fractal structure in general SOC models, this is not 
the case because the distribution of energy has to be of type like in the Zhang's model to get local contraction effects. 

We give now two simple examples of attractors which can be constructed "by hand" . 



2. One dimensional chain with E c 



1, N = 3,e = 0. 



For N — 3, each configuration X is a triplet {Xi, X3}. First note that only the mappings whose image intersect 
the cube [-^,-E c [ 3 are relevant for the asymptotic dynamics. Moreover, for E c < 1 each activation generates an 
avalanche, and the resulting configuration always contains a site with zero energy. This is an effect of projection onto 



the complementary set of the kernel of the product mapping, discussed in section II C 
The mappings (rather, their projection onto the faces of M) are respectively 
l 
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F 3 = 
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- 1 

! 


1 " 

I 


. 2 


4 . 


" 1 


" 





1 

2 . 


" 1 
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X 
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The Markov transition graph can be easily computed. Each legal transition occurs with probability ^ (activation 
of sites 1,2,3). To obtain the invariant set of the IFS, one must first notice that the three mappings F3,F5,F-j have 
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. These projection induce a 
tree structure for the invariant set. The maps send their domain of continuity onto the segments: 
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We can generate the invariant set by starting with the set : a U b U c. We show in Fig. 3 the initial branches a, b, c 
and their image under the five first iterates of the IFS. We have labeled the branches of the tree by their corresponding 
coding, for the three first iterates. One see, then how the tree structure is generated. 
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., >„ V.I' 1 



FIG. 1. Five first steps of iteration of the Iterated Function System , for E c = 1, N — 3,e — 0. The labeling of the branches 
of the tree, is a sequence, read from the left to the right, indicating the sequence of mappings applied to one of the initial 
branch a, b, c (right most symbol). The last picture is the plot of a trajectory on the attractor. 



3. One dimensional case where E c €]j^f , yr;]- 

This case is very atypical in the sense that the attractor is a finite set of points whose components have values • 
Indeed, let I be the set where each site has energy X = but at most one with energy j^. For E c e]y_f , yz^] 
the set I is the unique invariant set. Moreover, VB £ M, T n (B) — ► X, in the Haussdorf metric. Note that in this 
case we do not have invertibility, but each point has exactly N preimages. 

This behaviour is somehow pathological as it exists only in this range of E c value. For higher dimensions, we still 
do not know if there can be such an atomic invariant set. 

B. The probability distribution of avalanches size. 

In this part we derive a relation linking the sum of Lyapunov exponents and the probability of avalanche size. 
Then, we relate the fractal structure of the attractor to the critical exponent values. The basic ingredient is the 
Ledrappier- Young formula. However, as already said, the existence of a kernel in the standard Zhang's model makes 
the analysis somehow cumbersome. We therefore discuss first the non kernel case (e ^ 0) and comment only briefly 
on the modifications necessary to handle the limiting case e = 0. 

1. Average contraction rate. 

The key result is an exact formula linking the determinant of the basic maps to the total number s of overcritical 
sites in the corresponding avalanche. Namely we prove the following : 

Proposition 4 For £ K s one has 

detL k a = e s (21) 
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Proof We first show a propery about the relative distance of the overcritical sites in an avalanche. Let the avalanche 
be given by A = {Ak} 0<k<n where Ak is the set of overcritical sites at the kth step in the avalanche A. Denote 
by D(Ak) = {d(i,j) : i,j G Ak} the set of pairwise distances of the vertex set Ak- The proof is a straightforward 
consequence of the following lemma. 

Lemma 1 For any 7L d sublattice A, D(A k ) C 2. IN D(A k +i) C 2. IN 

Let denotes any path from i to j with no repetition of edges and I7I its length. From the general prop- 

erties of subsets of 7L d it follows that d(i,j) G 2.1SI =>• (7(2, j)| G 2. IN and, vice versa, if (7(2, j)| G 2. IN for some 
7 then d(i,j) E 2. IN. Let > 2 and i,j £ A k . We will show below, that provided e is sufficiently 

small, no site can be overcritical for two successive time steps. Assuming this for the moment, it follows that 
Ak+i C B(Ak, 1) = {v G A : d(v, Ak) = 1} since no site of Ak can be overcritical in the next step. Fix a path 7* by 
eliminating the first and last edge of 7 C Ak- 7* is a path between a vertex in B(i, 1) and B(j, 1), of length 7 — 2\. 
Since the pairwise distance in B(v, 1) are even for any vertex v E A it follows that for any pair of vertices from B(i, 1) 
to B(j, 1) there is an even length extension of 7* connecting those two vertices. This proves the lemma. 

We now show that, provided e is sufficiently small, a site cannot be overcritical in two successive time steps. For 
e = this is obvious. Assume now that e > 0. Let Ek be the maximal energy value of an overcritical site in the fcth 
step of an avalanche. For a given e we have to show that e.Ek < E c ,\/k. Clearly, Eq < E c + 1 . It is obvious that the 
maximal increase of energy on a site v can only happen if v has 2.d overcritical neighbours. In that case we have the 
following estimation: 



E k+ i < (l-e).E k + E c 

Iterating this expression we obtain: 



n_1 1 — (1 — f\ n 
E n <(l- e) n .E + E c . ]T(1 - eY = (1 - e) n .E + Ec. ^ '- 

i=o 6 
which has to be less than This holds provided : 

E c > (22) 

It follows from the lemma that two neighbours cannot be simultaneously overcritical during one avalanche provided 
e sufficiently small. One then gets the expression for the determinant by decomposing the matrix of the avalanche into 
one step matrices. The row corresponding to any overcritical site as only one non zero entry, the diagonal element e 
(nothing comes from the other overcritical sites at this time) while the columns corresponding to a non overcritical 
site has only one non zero entry, the diagonal element 1. Formula (21) follows. □ 

By using the ergodic theorem we get the log-average volume contraction which is also the sum of Lyapunov exponents 
as: 



1G 



Sn 

^ K = l °fJ € - ^2 sP N (s) = loge.s (23) 

i s—1 

where s is the average avalanche size and sn the maximal avalanche size. The formula relates the local volume con- 
traction to the average avalanche size. It connects therefore microscopic dynamical quantities (Lyapunov exponents) 
to a macroscopic observable (average avalanche size) . In particular it allows to establish a link between the Lyapunov 
spectrum and the critical exponents of the avalanche size distribution (see below and 0) 



2. Contraction versus expansion. 



The average contraction rate decreases with increasing E c . Indeed, the larger E c , the larger is the frequency of 
occurence of "trivial" avalanches where no relaxation occurs. They only display neutral directions in the phase space, 
with no contraction, and no contribution to the negative Lyapunov exponent. This can also be seen on formula ( |2l| ) 
: the larger E c , the smaller the average avalanche size. Therefore, for fixed N, there exists an E*(N) which is the 
unique E c value such that: 



loge.s + log(N) = 



(24) 



For E c < E*(N) the contraction dominates the expansion, while it is the opposite for E c > E*(N). Clearly, 
the invariant set structure is different in these two cases. On the one hand, for small E c values, the images of the 
domains V% are thin bands which are stretched slower than they contract. Therefore, they are expected not to overlap 
asymptotically and the invariant set has a Cantor structure with large gaps. On the other hand, when E c > E*(N), 
the successive images of the domains "P* fill more and more the phase space and the properties in conjecture 1 and 4 
should not hold. 

Note that, in this scheme, the Haussdorf dimension of A increases for increasing E c , E c < E*(N) and is likely to 
be constant when E c > E* (N) . 

The graph of E*(N) can easily be computed numerically. We give an example below, in a square lattice, for various 
values of e (Fig. |2|). Note that E*(N) increases with N. Therefore, one expects that, conjecture [l] || hold on larger 
and larger range of E c values, as N increases. 
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FIG. 2. E*(N) for various values of e. L is here the linear dimension of A (TV = L 2 ). 



3. Bounds on the critical exponent. 
In the invertible case the Ledrapier- Young formula implies: 

N 

^2\Xi\>logN (25) 

i=l 

Therefore, from formula ( p3| ) : 

This implies that the average avalanche size, s , has to diverge when N goes to infinity and that, in the thermody- 
namic limit, Pjv(s) tends to a distribution with an infinite mean- value. 

Furthermore, a reasonable assumption (supported by experiments) is that, for fixed N, the probability decreases 
with the avalanche size, namely: 

VN>N(e), P N (s)>P N (s + l) (27) 

In a certain way, this behaviour could be expected since the larger the avalanche, the more one has to impose 
conditions defining the corresponding domain of continuity, and the less the corresponding volume. However, this 
argument is not completely correct in general since one assumes some kind of absolute continuity of the invariant 
measure on the stable foliation (the probability of a domain decreases with its volume) . In particular it is completely 
false for Ec > 1 in the one dimensional chain, here the probability increases with s. 

Assuming that that there is indeed a power law, and that the system is invertible then one obtains: 

Pn(s) = K t < 2 (28) 

where /at(s) is a cut-off function accounting for finite size effects 

Therefore, eq. (26) gives the scaling of the power law and bounds for the critical exponent r. In particular, if we 
assume that Pn{s) converges to some limit P*(s) as N — > +oo, then P*(s) — -pr, t G [1, 2]. 

4- The value of r and the fractality of the support of the invariant measure. 

The Ledrappier- Young formula gives a direct way to check the "fractality" of the support of the invariant measure 
in the invertible case. For 



8 In the SOC litterature, the Finite Scaling Assumption leads to write the probability distribution of avalanche size Pjv(s) as 
Pn(s) = s~ T Q(s.L~ /3 ) where Q accounts for cut-off effects. The exponents (t, (3) are believed to characterize the universality 
class of the model. Note that r, j3 depend a priori from Ec. 
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\Y^Xr\>logN (29) 

i 

if some partial dimensions <7, are not integers (< mi). 

Suppose that the measure [i is absolutely continuous on the stable foliation. Then from eq. ( |i~7| ) we get equality 



in (26), implying that the sum diverges logarithmically with N. Furthermore, the maximal avalanche size scales like 
sn w N 13 / 2 , implying that logsN ~ logN . Then s diverges logarithmically with sn suggesting a critical exponent 
t = 2. More generally, we get the same result if the fractal set is homogeneous (in sense that all partial Haussdorf 
dimensions are equal or can be bounded from below as N — ► oo). 

However, one does not expect the fractal to be homogeneous. It is indeed clear that the contraction is not uniform 
in the phase space. For e = the kernel directions produce infinite contraction. In figure [l] they are the directions 
transverse to the "branches" of the attractor, which project the dynamics on the tree, in one time step. On the 
other hand, the directions "parallel" to the branches produce finite contraction. As a corrolary, the partial Hausdorff 
dimensions of the invariant set are zero transversally to the attractor while they are finite along the branches. When 
e is small, there are still directions producing high contractions, those which give the kernel directions as e — ► 0. This 
effect is reflected in the Lyapunov spectrum where one detects two parts in the spectrum (see Figure ^|). 
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FIG. 3. Spectrum of negative Lyapunov exponents for N = 49, Ec = 3, e = 0.1. 



Therefore, the sum (resp. s) must diverge faster than logarithmically with N and strict inequality holds. This 
implies that the critical exponent r < 2. The measured exponent is indeed always strictly lower than 2 pi] . Note 
that under the finite size scaling hypothesis s behaves as N^-^i and that r is indeed lower than 2 iff strict inequality 
holds in @. 

An explicit formula linking the Hausdorff dimension and the critical exponents (r, (3) can be obtained through the 
Ledrappier- Young formula. This will be treated in a separated paper H . 



19 



5. The case e = 



One would like to obtain an equality like (|2g) also in this case. However, setting e — leads to a zero determinant an 
hence infinite Lyapunov exponents (more precisely the Lyapunov exponents are not defined on the whole configuration 
space). But restricted to Z S (X) one can still compute finite Lyapunov exponents and the sum is just the determinant 
of the matrix restricted to the stable space E S (X) at some point X in the domain of the map. 

If one has a furhter inequality like: 

det{T k a \ £s{±) )>C s (30) 

VK G K s and < C < 1 one still gets the same kind of estimates for the expected avalanche size like in the e ^ 
case. The details are quite cumbersome and will be given in a forthcoming paper. 

C. Phase transitions. 

The domains of continuity hA\ are bounded by hyperplanes, which are moving when E c varies. In general, a small 
variation in E c does not lead to structural changes in the dynamics, if all these hyperplanes are intersecting the interior 
of Ai. In this case, the structure of the transition graph is not modified. Moreover, the corresponding mapping T k 
does not change under this motion. More precisely, changes in E c just change the shape of M.\ but not the matrix of 
the mapping T*. 

However, for some E c values, some hyperplanes have intersection only with d A4. This implies that a small change 
in E c can push these hyperplanes outside M. Hence the corresponding transition graph changes in structure. As far 
as the asymptotic dynamics and therefore, the invariant distribution is dependent on the graph structure, we expect 
changes in the SOC picture when crossing these critical E c values. This effect has already been reported elsewhere for 
the one dimensional Zhang's model and arises also in two dimensions where Pn(s) is not a power law for E c << 1 

§• 

In fact, one can easily figure out that at least the limiting cases E c — > oo and E c — > are completely different. 
For E c — > oo relaxation events are more and more seldom. One obtains kind of a frozen state where energy increases 
(on average) monotonously with some rare (but large) avalanches. Moreover, the asymptotic energy distribution is 
sensitive to the initial conditions (loss of ergodicity). Furthermore, the attractor as a large Haussdorf dimension. 

On the other hand, for E c — > 0, each activation generates a very large avalanche (that has to reflect many times 
on the boundary before it has lost enough energy to stop). This implies larger and larger contraction, and therefore 
the sum of Lyapunov exponents decreases to — oo. As a corollary of Ledrappier- Young formula the partial fractal 
dimensions have to go to zero in order to maintain the product equal to logN. 

IV. CONCLUSION. 

We have shown that certain classes of models of SOC like the Zhang's model fit naturally into a well known class 
of dynamical systems. Especially for the question of asymptotic energy distribution, observables distribution, ergod- 
icity, this seems to be a proper point of view. Furthermore it seems likely to exhibit close relationship between the 
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probability of the size of avalanches and the fractality of the attractor. 

There are many questions for further investigations. We list a few of them. 

1. Development of a thermodynamic formalism and it's linkage to the SOC quantities. It should 
be possible to extrapolate this formalism to the case of arbitrary (hyperbolic) SOC-system. Moreover, phase 
transitions should correspond to changes in the invariant measure of maximal entropy (loss of analyticity of 
the topological presssure). One expects that in a proper formulation E c should play the role of an inverse 
temperature. 

2. Dimension spectrum of the attractor and Lyapunov spectrum. We have outlined above the crucial 
role played by Lyapunov exponents (accounting for energy transport) and the link one can establish with the 
equilibrium state and the critical exponents. The full developments of this point will be published elsewhere j^] 

3. Nonuniform distribution rules. As outlined in the paper most of our results carry on if one does not choose 
a uniform activation measure, because one still has a good measure as an equilibrium state. On the other hand, 
activating with a degenerate probability distribution (for example activating always the same site) will lead 
to different results. Activating sites periodically with different period will allow to sample the periodic orbits 
structure of the global attractor, which are dense. 

4. Thermodynamic limit for fixed E c and N — > oo and the limit E c — > oo (N fixed). In these both cases 
one loses the hyperbolic structure. 

5. Smooth thresholds. Some modification of Zhang's model have been proposed, in particular to treat this 
model in the continuum limit by an anomalous diffusion equation |l^,^|]. In this case the Heaviside function 
corresponding to the sharp threshold at E c is smoothed out by some continuous function. The nice effect of this 
change in our description is that it removes the singularity set. On the other hand, the system is expected to 
have still a nice hyperbolic structure (though non uniform) where smooth local stable manifolds exist for almost 
all points. Pesin theory |2^] should apply in this context. 

6. The case where 5X is random. In the usual Zhang's model, the energy activation quantum SX is not 
a constant but is a random variable. This situation can be treated in the framework of random hyperbolic 
dynamical system. 

As a conclusion we would like to outline that the study of SOC-models with tools from dynamical system theory 
will certainly not solve all questions in this context. In the belief of the authors it is mainly useful for the study 
of fairly general structure properties of the models. It is also clear that the complexity of the underlying transition 
graph on which the model is defined will become of crucial importance for some questions. 
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List of Figures 



1 Five first steps of iteration of the Iterated Function System , for E c = 1, N = 3, e = 0. The labeling 
of the branches of the tree, is a sequence, read from the left to the right, indicating the sequence of 
mappings applied to one of the initial branch a,b,c (right most symbol). The last picture is the plot of 
a trajectory on the attractor 

2 E*(N) for various values of e. L is here the linear dimension of A (N = L 2 ) 

3 Spectrum of negative Lyapunov exponents for N = 49, Ec — 3, e = 0.1 
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